#######################################
######Masculinity Threat Analysis######
#######################################
###############10/10/25##################
#######################################

####Analysis Script for Gothreau & Haas, forthcoming in Journal of Experimental Political Science###
####A Replication and Extension of Willer et al. (2013) Overdoing Gender: A Test of the Masculine Overcompensation Thesis###

library(haven)
library(ggplot2)
library(dplyr)
library(broom)
library(officer)
library(flextable)
library(forcats)
library(psych)
library(survey)
library(readr)
library(tidyr)

#Set Working Directory#
MTData <- read.csv("CleanMTData.csv")

table(MTData$P_EXP, MTData$gender_binary)

#Create datasets by gender
Data_female <- subset(MTData, gender_binary == "Women")
Data_male <- subset(MTData, gender_binary == "Men")

#Set gender and experimental conditions as factors 
MTData$P_EXP <- factor(MTData$P_EXP, levels=c("1", "2", "3", "4"), labels = c("threat", "no threat", "constrained", "PCKT"))
MTData$gender_binary <- factor(MTData$Q_CDGDRDESC_17, levels = c("1", "2"), labels = c("Men", "Women"))

#WEIGHTED T-TESTS 
#List of dependent variables
dvs <- c("IRAQ_WAR_rescaled", "GAY_RIGHT_rescaled", "SUV_DESIRABILITY_rescaled", "SUV_PAY_rescaled")

run_weighted_ttests <- function(data, label) {
  experimental_groups <- setdiff(unique(data$P_EXP), "no threat")
  results <- list()
  
  for (dv in dvs) {
    if (!dv %in% names(data) || !is.numeric(data[[dv]])) {
      warning(paste("Skipping", dv, "- missing or not numeric"))
      next
    }
    
    for (group in experimental_groups) {
      # Subset data to include only "no threat" and one comparison group
      subset_data <- data %>%
        filter(P_EXP %in% c("no threat", group)) %>%
        select(P_EXP, WEIGHT, all_of(dv)) %>%
        na.omit()
      
      if (nrow(subset_data) < 3) next  # skip if too few cases
      
      # Define survey design
      design <- svydesign(ids = ~1, weights = ~WEIGHT, data = subset_data)
      
      # Run weighted t-test
      formula <- as.formula(paste(dv, "~ P_EXP"))
      ttest <- svyttest(formula, design)
      
      # Extract test statistics
      t_val <- ttest$statistic
      p_val <- ttest$p.value
      
      # Get weighted means
      group_means <- svyby(as.formula(paste("~", dv)), ~P_EXP, design, svymean, na.rm = TRUE)
      mean_no_threat <- group_means[group_means$P_EXP == "no threat", 2]
      mean_treat <- group_means[group_means$P_EXP == group, 2]
      
      result <- data.frame(
        gender = label,
        DV = dv,
        experimental_group = group,
        mean_no_threat = round(mean_no_threat, 2),
        mean_treat = round(mean_treat, 2),
        t = round(t_val, 2),
        p = p_val
      )
      
      result$sig <- ifelse(result$p < .001, "***",
                           ifelse(result$p < .01, "**",
                                  ifelse(result$p < .05, "*", "")))
      
      results[[paste(label, dv, group, sep = "_vs_")]] <- result
    }
  }
  
  bind_rows(results)
}

ttests_female <- run_weighted_ttests(Data_female, "female")
ttests_male   <- run_weighted_ttests(Data_male, "male")

all_ttests <- bind_rows(ttests_female, ttests_male)

######Replicate Willer et al. Study 1 Table (Table 1 in Main Text)

# Keep only comparisons where experimental group is "threat"
threat_table <- all_ttests %>%
  filter(experimental_group == "threat") %>%
  mutate(
    comparison = paste0("t = ", t, sig, ", p = ", sprintf("%.3f", p)),
    DV = factor(DV, levels = dvs)
  ) %>%
  select(gender, DV, mean_treat, mean_no_threat, comparison)

ft <- flextable(threat_table)
ft <- footnote(
  ft,
  i = 1, j = 1,
  value = as_paragraph("* p < .05, ** p < .01, *** p < .001"),
  ref_symbols = "†",
  part = "body"
)
ft <- add_footer_lines(ft, "* p < .05, ** p < .01, *** p < .001; † indicates significance in original study")
ft <- autofit(ft)

doc <- read_docx()
doc <- body_add_par(doc, "T-Test Results: Threat vs. No Threat Condition", style = "heading 1")
doc <- body_add_flextable(doc, ft)
print(doc, target = "Study1_replication.docx")

######Replicate Table C.1 in the Appendix 
#################P-Value Adjustment for Multiple Comparison##############

dvs_adjusted <- c(
  "IRAQ_WAR_rescaled", "GAY_RIGHT_rescaled", "SUV_DESIRABILITY_rescaled", "SUV_PAY_rescaled",
  "SYSTEM_JUSTIFICATION_rescaled", "DOMINANCE_scale", "TRADITIONALISM_scale",
  "CONSERVATISM_rescaled", "IMMIGRATION_rescaled", "MARIJUANA_rescaled",
  "ELECTRIC_CAR_rescaled", "TRANS_RIGHT_rescaled", "HIRINGPOLICY_WOMEN_rescaled"
)

run_ttests_bh <- function(data, gender_label) {
  # comparisons that actually exist besides "no threat"
  px <- data$P_EXP
  comparisons <- setdiff(unique(as.character(px)), "no threat")
  comparisons <- comparisons[!is.na(comparisons) & nzchar(comparisons)]
  
  results <- list()
  
  for (dv in dvs_adjusted) {
    for (cond in comparisons) {
      sub_data <- data %>%
        filter(P_EXP %in% c("no threat", cond)) %>%
        select(P_EXP, WEIGHT, all_of(dv)) %>%
        drop_na()
      
      # need at least 3 obs and both groups present
      if (nrow(sub_data) < 3) next
      
      # force two-level factor in the right order
      sub_data <- sub_data %>%
        mutate(P_EXP = factor(P_EXP, levels = c("no threat", cond))) %>%
        filter(!is.na(P_EXP))
      
      if (nlevels(sub_data$P_EXP) != 2) next
      # ensure both groups have at least 1 obs
      if (any(table(sub_data$P_EXP) == 0)) next
      
      # survey design
      design <- svydesign(ids = ~1, weights = ~WEIGHT, data = sub_data)
      fmla <- as.formula(paste(dv, "~ P_EXP"))
      
      # run svyttest; skip on error
      test <- try(svyttest(fmla, design), silent = TRUE)
      if (inherits(test, "try-error")) next
      
      # group means
      means <- svyby(as.formula(paste0("~", dv)), ~P_EXP, design, svymean, na.rm = TRUE)
      mean_no_threat <- means %>% filter(P_EXP == "no threat") %>% pull(2)
      mean_treat <- means %>% filter(P_EXP == cond) %>% pull(2)
      
      results[[paste(gender_label, dv, cond, sep = "_")]] <- tibble(
        gender = gender_label,
        DV = dv,
        experimental_group = cond,
        mean_no_threat = round(mean_no_threat, 3),
        mean_treat = round(mean_treat, 3),
        t = round(unname(test$statistic), 2),
        p = unname(test$p.value)
      )
    }
  }
  
  if (length(results) == 0) {
    # return an empty tibble with expected columns so downstream code won’t fail
    return(tibble(
      gender = character(),
      DV = character(),
      experimental_group = character(),
      mean_no_threat = numeric(),
      mean_treat = numeric(),
      t = numeric(),
      p = numeric(),
      p_bh = numeric(),
      stars = character(),
      comparison = character()
    ))
  }
  
  final <- bind_rows(results)
  
  # Adjust p-values using Benjamini-Hochberg (BH) within gender
  final <- final %>%
    group_by(gender) %>%                        # adjust within gender
    mutate(p_bh = p.adjust(p, method = "BH")) %>%
    ungroup() %>%
    mutate(
      stars = case_when(
        p_bh < .001 ~ "***",
        p_bh < .01  ~ "**",
        p_bh < .05  ~ "*",
        TRUE ~ ""
      ),
      comparison = paste0("t = ", t, stars, ", p = ", sprintf("%.3f", p_bh))
    )
  
  final
}

# Run
ttests_female <- run_ttests_bh(Data_female, "Female")
ttests_male   <- run_ttests_bh(Data_male, "Male")
all_results   <- bind_rows(ttests_female, ttests_male)

# If empty, you can bail out early or proceed to a blank table
if (nrow(all_results) == 0) {
  warning("No valid t-tests were computed. Check that P_EXP contains 'no threat' and at least one other condition with data for each DV.")
}

table_data <- all_results %>%
  select(gender, DV, experimental_group, mean_no_threat, mean_treat, comparison)

ft <- flextable(table_data) %>%
  set_header_labels(
    gender = "Gender",
    DV = "Outcome",
    experimental_group = "Condition",
    mean_no_threat = "Mean (No Threat)",
    mean_treat = "Mean (Treatment)",
    comparison = "Comparison (BH Adjusted)"
  ) %>%
  width(j = "gender", width = 1.0) %>%
  width(j = "DV", width = 1.2) %>%
  width(j = "experimental_group", width = 1.5) %>%
  width(j = "mean_no_threat", width = 1.5) %>%
  width(j = "mean_treat", width = 1.5) %>%
  width(j = "comparison", width = 2.3) %>%
  # set_table_properties: prefer width = 1 (inches) OR layout="autofit", not both fighting each other
  set_table_properties(layout = "autofit") %>%
  align(j = "comparison", align = "left", part = "all") %>%
  autofit() %>%
  add_footer_lines("* p < .05, ** p < .01, *** p < .001; p-values Benjamini–Hochberg adjusted")

doc <- read_docx()
doc <- body_add_par(doc, "Weighted T-Test Results by Gender and Condition", style = "heading 1")
doc <- body_add_flextable(doc, ft)
print(doc, target = "TTest_BH_Adjusted.docx")

###Replication of Willer et al Study 2 
#List of dependent variables
dvs <- c("SYSTEM_JUSTIFICATION_rescaled", "DOMINANCE_scale", "TRADITIONALISM_scale", "CONSERVATISM_rescaled")

run_weighted_ttests <- function(data, label) {
  experimental_groups <- setdiff(unique(data$P_EXP), "no threat")
  results <- list()
  
  for (dv in dvs) {
    if (!dv %in% names(data) || !is.numeric(data[[dv]])) {
      warning(paste("Skipping", dv, "- missing or not numeric"))
      next
    }
    
    for (group in experimental_groups) {
      # Subset data to include only "no threat" and one comparison group
      subset_data <- data %>%
        filter(P_EXP %in% c("no threat", group)) %>%
        select(P_EXP, WEIGHT, all_of(dv)) %>%
        na.omit()
      
      if (nrow(subset_data) < 3) next  # skip if too few cases
      
      # Define survey design
      design <- svydesign(ids = ~1, weights = ~WEIGHT, data = subset_data)
      
      # Run weighted t-test
      formula <- as.formula(paste(dv, "~ P_EXP"))
      ttest <- svyttest(formula, design)
      
      # Extract test statistics
      t_val <- ttest$statistic
      p_val <- ttest$p.value
      
      # Get weighted means
      group_means <- svyby(as.formula(paste("~", dv)), ~P_EXP, design, svymean, na.rm = TRUE)
      mean_no_threat <- group_means[group_means$P_EXP == "no threat", 2]
      mean_treat <- group_means[group_means$P_EXP == group, 2]
      
      result <- data.frame(
        gender = label,
        DV = dv,
        experimental_group = group,
        mean_no_threat = round(mean_no_threat, 2),
        mean_treat = round(mean_treat, 2),
        t = round(t_val, 2),
        p = p_val
      )
      
      result$sig <- ifelse(result$p < .001, "***",
                           ifelse(result$p < .01, "**",
                                  ifelse(result$p < .05, "*", "")))
      
      results[[paste(label, dv, group, sep = "_vs_")]] <- result
    }
  }
  
  bind_rows(results)
}

ttests_female <- run_weighted_ttests(Data_female, "female")
ttests_male   <- run_weighted_ttests(Data_male, "male")

all_ttests <- bind_rows(ttests_female, ttests_male)


######Replicate Willer et al. Study 2 Table (Table 2 in Main Text)

# Keep only comparisons where experimental group is "threat"
threat_table <- all_ttests %>%
  filter(experimental_group == "threat") %>%
  mutate(
    comparison = paste0("t = ", t, sig, ", p = ", sprintf("%.3f", p)),
    DV = factor(DV, levels = dvs)
  ) %>%
  select(gender, DV, mean_treat, mean_no_threat, comparison)

ft <- flextable(threat_table)
ft <- footnote(
  ft,
  i = 1, j = 1,
  value = as_paragraph("* p < .05, ** p < .01, *** p < .001"),
  ref_symbols = "†",
  part = "body"
)
ft <- add_footer_lines(ft, "* p < .05, ** p < .01, *** p < .001; † indicates significance in original study")
ft <- autofit(ft)

doc <- read_docx()
doc <- body_add_par(doc, "T-Test Results: Threat vs. No Threat Condition", style = "heading 1")
doc <- body_add_flextable(doc, ft)
print(doc, target = "Study2_replication.docx")





















































